% Function to estimate the drag of an aircraft
% with quadratic (parabolic) drag polar
% Can be used stand-alone or with various performance programs

% Programmer: Toni Trani 
% creation date: 09/10/03
% modification 1 = added input variable "phi" to calculate drag for a given
% bank angle (09/12/04)

% Outputs:

% drag                  = total aircraft drag (N)
% LD                    = lift to drag ratio (dimensionless)
% lift_coefficient    = lift coefficient (dimensionless)
% drag_coefficient = total drag coefficient (dimensionless)

% Inputs:

% mach                = true mach number (dimensionless)
% rho                   = air density (kg/cu. meter)
% a                      = speed of sound (m/s-s)
% mass                 = aircraft mass (kg)
% macht               = table of mach numbers for Cdo computation (dim)
% Cdoct               = table of values of zero lift drag coefficient - Cdo (dim)
% phi                   = bank angle (radians)

function [drag,LD,lift_coefficient,drag_coefficient]=drag_03(mach,rho,a,mass,macht,Cdoct,phi)

global A e S g          % aspect ratio, wing efficiency factor, wing area and gravity constant

% start computations
	rhos = 1.225;                             % value of density at sea level conditions (kg/cu.m.)    
	vtas = mach .* a;                        % true airspeed (m/s)
	q = 0.5 * rho .* vtas.^2;              % dynamic pressure (N/sq. meter)

	% Estimates the cruise CL (lift coefficient) from basic lift equation

	lift_coefficient = g*mass ./ (q .* S .* cos(phi));

	% Estimates drag coefficient (first interpolate to get Cdo, then find
	% total drag coefficient)

	cdo =  	interp1(macht,Cdoct,mach,'pchip');                                                    % zero lift drag coefficient (dim)
	drag_coefficient 	= 	cdo + lift_coefficient.^2 ./ (pi.* A.*e);                 % total drag coefficient	(dim)								

	% Estimate the total vehicle drag
	
	LD = lift_coefficient ./ drag_coefficient ;                                              % Lift to drag ratio (dim)
	drag =  S.*q.* drag_coefficient ;                                                          % total drag (N)
